library(tidyverse)
library(lmerTest)pal_process = wesanderson::wes_palette("Zissou1", 6, "continuous")source("load_data.R")# standardize and winsorize
betas_std = betas %>%
group_by(roi, session) %>%
mutate(meanPE_std = scale(meanPE, center = TRUE, scale = TRUE),
meanPE_std = ifelse(meanPE_std > 3, 3,
ifelse(meanPE_std < -3, -3, meanPE_std))) %>%
ungroup()
dots_std = dots %>%
group_by(map, test, mask, session) %>%
mutate(dotProduct_std = scale(dotProduct, center = TRUE, scale = TRUE),
dotProduct_std = ifelse(dotProduct_std > 3, 3,
ifelse(dotProduct_std < -3, -3, dotProduct_std))) %>%
ungroup()
# join data frames
dataset = full_join(betas_std, dots_std, by = c("subjectID", "con", "process", "condition", "control", "session")) %>%
mutate(subjectID = as.character(subjectID)) %>%
ungroup() %>%
mutate(subjectID = as.factor(subjectID),
condition = as.factor(condition),
control = as.factor(control),
roi = as.factor(roi),
process = as.factor(process),
test = as.factor(test))ema_tidy = ema %>%
select(subjectID, desire_pres0.prop, desire_str.m, goal_conf.m, resist.m, enact_prop, amount_Recoded.m, hunger.m) %>%
unique()
model_data = dataset %>%
filter(((!grepl("neuralsig", map) & test == "association" & mask == "masked") | (grepl("neuralsig", map) & mask == "unmasked")) &
session == "all" & control == "nature" & condition == "snack") %>%
unique() %>%
group_by(process, subjectID) %>%
mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE),
processPE = paste0(process, "_ROI"),
processPEV = paste0(process, "_PEV")) %>%
ungroup() %>%
select(-c(xyz, meanPE, meanPE_std, sdPE, dotProduct, map, process)) %>%
spread(processPEV, dotProduct_std) %>%
spread(processPE, meanProcessPEstd) %>%
group_by(subjectID) %>%
fill(everything(), .direction = "down") %>%
fill(everything(), .direction = "up") %>%
select(-c(roi, craving_ROI, craving_regulation_ROI, test, mask)) %>%
unique() %>%
mutate(balance_ROI = cognitive_control_ROI - reward_ROI) %>%
right_join(., ema_tidy, by = "subjectID") %>%
na.omit()model_data %>%
select(subjectID, contains("ROI"), desire_pres0.prop, desire_str.m, goal_conf.m, resist.m, enact_prop, amount_Recoded.m) %>%
unique() %>%
gather(neuro_var, neuro_val, contains("ROI")) %>%
mutate(neuro_var = factor(neuro_var, levels = c("balance_ROI", "cognitive_control_ROI",
"reward_ROI", "value_ROI"))) %>%
gather(ema_var, ema_val, contains("m"), contains("prop")) %>%
group_by(ema_var) %>%
do({
plot = ggplot(., aes(neuro_val, ema_val, color = neuro_var)) +
geom_point(alpha = .5) +
geom_smooth(method = "lm", alpha = .1, fullrange = TRUE) +
scale_color_manual(name = "", values = pal_process) +
labs(x = "\n mean parameter estimate", y = paste0(.$ema_var[[1]], "\n")) +
theme_minimal(base_size = 14)
print(plot)
data.frame()
})model_data %>%
select(subjectID, contains("PEV"), desire_pres0.prop, desire_str.m, goal_conf.m, resist.m, enact_prop, amount_Recoded.m) %>%
unique() %>%
gather(neuro_var, neuro_val, contains("PEV")) %>%
mutate(neuro_var = factor(neuro_var, levels = c("cognitive_control_PEV", "reward_PEV", "value_PEV",
"craving_PEV", "craving_regulation_PEV"))) %>%
gather(ema_var, ema_val, contains("m"), contains("prop")) %>%
group_by(ema_var) %>%
do({
plot = ggplot(., aes(neuro_val, ema_val, color = neuro_var)) +
geom_point(alpha = .5) +
geom_smooth(method = "lm", alpha = .1, fullrange = TRUE) +
scale_color_manual(name = "", values = pal_process[2:6]) +
labs(x = "\n mean parameter estimate", y = paste0(.$ema_var[[1]], "\n")) +
theme_minimal(base_size = 14)
print(plot)
data.frame()
})options(na.action = "na.fail")
data_ctrl = caret::trainControl(method = "repeatedcv", number = 5, repeats = 5)desire_prop_roi_full = lm(desire_pres0.prop ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data)
(desire_prop_roi_mods = MuMIn::dredge(desire_prop_roi_full))MuMIn::get.models(desire_prop_roi_mods, subset = TRUE) %>%
tibble() %>%
rename("model" = ".") %>%
mutate(tidied = purrr::map(model, broom::tidy),
model_num = row_number()) %>%
select(model_num, tidied) %>%
unnest()best_desire_prop_roi = caret::train(desire_pres0.prop ~ reward_ROI,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)desire_prop_mam_full = lm(desire_pres0.prop ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data)
(desire_prop_mam_mods = MuMIn::dredge(desire_prop_mam_full))MuMIn::get.models(desire_prop_mam_mods, subset = TRUE) %>%
tibble() %>%
rename("model" = ".") %>%
mutate(tidied = purrr::map(model, broom::tidy),
model_num = row_number()) %>%
select(model_num, tidied) %>%
unnest()best_desire_prop_mam = caret::train(desire_pres0.prop ~ value_PEV,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)desire_prop_null = lm(desire_pres0.prop ~ 1, data = model_data)
desire_prop_combined = caret::train(desire_pres0.prop ~ reward_ROI + value_PEV,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)
AIC(desire_prop_null, best_desire_prop_roi$finalModel, best_desire_prop_mam$finalModel, desire_prop_combined$finalModel) %>%
rownames_to_column() %>%
arrange(AIC)Summarize the overall best fitting
summary(best_desire_prop_mam$finalModel)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.27760 -0.12561 -0.01319 0.10406 0.37274
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.33844 0.01682 20.119 <2e-16 ***
## value_PEV 0.02656 0.01382 1.922 0.0576 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1494 on 95 degrees of freedom
## Multiple R-squared: 0.03743, Adjusted R-squared: 0.0273
## F-statistic: 3.694 on 1 and 95 DF, p-value: 0.05761
best_desire_prop_mam## Linear Regression
##
## 97 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times)
## Summary of sample sizes: 77, 78, 78, 78, 77, 78, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.1496044 0.09099185 0.1255609
##
## Tuning parameter 'intercept' was held constant at a value of TRUE
sd(best_desire_prop_mam$resample$Rsquared)## [1] 0.102751
sd(best_desire_prop_mam$resample$RMSE)## [1] 0.01693802
enact_prop_roi_full = lm(enact_prop ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data)
(enact_prop_roi_mods = MuMIn::dredge(enact_prop_roi_full))MuMIn::get.models(enact_prop_roi_mods, subset = TRUE) %>%
tibble() %>%
rename("model" = ".") %>%
mutate(tidied = purrr::map(model, broom::tidy),
model_num = row_number()) %>%
select(model_num, tidied) %>%
unnest()best_enact_prop_roi = caret::train(enact_prop ~ reward_ROI,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)enact_prop_mam_full = lm(enact_prop ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data)
(enact_prop_mam_mods = MuMIn::dredge(enact_prop_mam_full))MuMIn::get.models(enact_prop_mam_mods, subset = TRUE) %>%
tibble() %>%
rename("model" = ".") %>%
mutate(tidied = purrr::map(model, broom::tidy),
model_num = row_number()) %>%
select(model_num, tidied) %>%
unnest()best_enact_prop_mam = caret::train(enact_prop ~ cognitive_control_PEV + craving_regulation_PEV + reward_PEV + value_PEV,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)enact_prop_null = lm(enact_prop ~ 1, data = model_data)
enact_prop_combined = caret::train(enact_prop ~ reward_ROI + cognitive_control_PEV + craving_regulation_PEV + reward_PEV + value_PEV,
data = model_data,
trControl = data_ctrl,
method = "lm",
na.action = na.omit)
AIC(enact_prop_null, best_enact_prop_roi$finalModel, best_enact_prop_mam$finalModel, enact_prop_combined$finalModel) %>%
rownames_to_column() %>%
arrange(AIC)Summarize the overall best fitting
summary(best_enact_prop_mam$finalModel)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.61127 -0.13058 0.00977 0.10827 0.49466
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.47037 0.02893 16.261 < 2e-16 ***
## cognitive_control_PEV -0.06868 0.02749 -2.498 0.01426 *
## craving_regulation_PEV 0.03779 0.02324 1.626 0.10731
## reward_PEV 0.19511 0.06939 2.812 0.00602 **
## value_PEV -0.12959 0.06221 -2.083 0.04001 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2121 on 92 degrees of freedom
## Multiple R-squared: 0.1349, Adjusted R-squared: 0.09732
## F-statistic: 3.587 on 4 and 92 DF, p-value: 0.009164
best_enact_prop_mam## Linear Regression
##
## 97 samples
## 4 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times)
## Summary of sample sizes: 78, 77, 78, 77, 78, 78, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.214186 0.1087687 0.1628077
##
## Tuning parameter 'intercept' was held constant at a value of TRUE
sd(best_enact_prop_mam$resample$Rsquared)## [1] 0.102848
sd(best_enact_prop_mam$resample$RMSE)## [1] 0.03547659